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Precise measurements of the anisotropies in the cosmic microwave background enable us to do 
an accurate study on the form of the primordial power spectrum for a given set of cosmological 
parameters. In a previous paper [l|, we implemented an improved (error sensitive) Richardson-Lucy 
deconvolution algorithm on the measured angular power spectrum from the first year of WMAP 
data to determine the primordial power spectrum assuming a concordance cosmological model. 
This recovered spectrum has a likelihood far better than a scale invariant, or, 'best fit' scale free 
spectra {A\n£. ~ 25 w.r.t. Harrison Zeldovich, and, Aln£ « 11 w.r.t. power law with Ua = 0.95). 
In this paper we use Discrete Wavelet Transform (DWT) to decompose the local features of the 
recovered spectrum individually to study their effect and significance on the recovered angular 
power spectrum and hence the likelihood. We show that besides the infra-red cut off at the horizon 
scale, the associated features of the primordial power spectrum around the horizon have a significant 
effect on improving the likelihood. The strong features are localized at the horizon scale. 



Observational comparison of cosmological models based on structure formation in the universe necessarily depends 
on the assumed initial conditions describing the primordial seed perturbations. Inflation, besides resolving a number 
of problems of classical Big Bang theory, provides us a mechanism for generating these correlated primordial per- 
turbations 0, S 0] ■ Precision measurements of anisotropies in the cosmic microwave background, and also of the 
clustering of large scale structure, suggest that the primordial density perturbation is dominantly adiabatic and has 
a nearly scale invariant spectrum This is in good agreement with most simple inflationary scenarios which 

predict power law or scale invariant forms of the primordial perturbation. The data have also been used widely to 
put constraints on different parametric forms of primordial spectrum, mostly motivated by inflation 0, d, [l0| . 

However, despite the strong theoretical appeal and simplicity of a featureless primordial spectrum, it is important to 
determine the shape of the primordial power spectrum directly from observations with minimal theoretical bias. Many 
model independent searches have been made to look for features in the CMB primordial power spectrum [ll|, [l^, [isl . 



deconvolved directly from the angular power spectrum of CMB anisotropy measured by WMAP using an improved 
implementation of the Richardson-Lucy algorithm The most prominent feature of the recovered spectrum is 
a sharp infra-red cut-off on the horizon scale. It also has a localized excess just above the cut-off which leads to 
great improvement of likelihood over the simple monotonic forms of model infra-red cut-off spectra considered widely 
in the post-WMAP-1 literature. Interestingly, similar features were also detected by method of regularized least 
squares [l9| . The signiflcant improvement in the likelihood clearly shows the importance of features in the primordial 
power spectrum. 

The goal of the present paper is to demonstrate the application of wavelets in identifying the statistically significant 
features in a deconvolved power spectrum. We use Discrete Wavelet Transform (henceforth DWT) to identify features 
in the recovered primordial power spectrum at different resolutions and at different locations in k space. Starting 
from the coarsest primordial power spectrum we systematically add variations on different resolutions and obtain 
the angular power spectrum. We then compute the likelihood of the reconstructed primordial power spectrum by 
comparing the angular power spectrum with the WMAP data. The improvement in the likelihood allows us to quantify 
the significance of different features. 



The Richardson-Lucy (RL) algorithm was developed and is widely used in the context of image reconstruction in 
astronomy [2^,[3^. However, the method has also been successfully used in cosmology, to deproject the 3-D correlation 



I. 



INTRODUCTION 




II. RICHARDSON-LUCY DECONVOLUTION METHOD 



2 



function and power spectrum from the measured 2-D angular correlation and 2-D power spectrum [3l|, [S^l ■ 

The angular power spectrum, C;, is a convolution of the initial power spectrum P{k) generated in the early universe 
with a radiative transport kernel, G{l,k), that is determined by the values of the cosmological parameters. In our 
application, we solve the inverse problem of determining the primordial power spectrum, -P(fc), from the measured 
angular power spectrum, C; , using the relation 



G =^Ga,fc,)-P(fc^)- (1) 

i 

In the above equation, the 'target' measured angular power spectrum, Ci = Cf , is the data given by observations, 
and the radiative transport kernel. 

Ah. 

G{l,h)^—^\ATi{h,fjo)\\ (2) 

encodes the response of the present multipoles of the CMB perturbed photon distribution function ATi{ki,riQ) to unit 
of power per logarithm interval of wavenumber, k, in the primordial perturbation spectrum. The kernel G{1, k) is 
completely fixed by the cosmological parameters of the 'base ' cosmological model. Obtaining P{k) from the measured 
C/, for a given G{l,k), is clearly a deconvolution problem. An important feature of the problem is that Gf, G{l,k) 
and P{k) are all positive definite. 

In Ref. [l| we (A.S. and T.S.) employ an improved RL method to solve the inverse problem for P{k) in Eq. ([T|). 
The advantage of RL method is that positivity of the recovered P{k) is automatically ensured, given G{l,k) is 
positive definite and C/'s are positive. The RL method is readily derived from elementary probability theory on 
distributions [2^ and is an iterative method which can be neatly encoded into a simple recurrence relation. The 
power spectrum p(*+^)(A:) recovered after iteration {i + 1) is given by 



p(^+i)(fc)-pW(fc) ^pW(fc)^G(/,fc)^^^^^ (3) 

I 

where Gf is the measured data (target) and G\^^ is the angular power spectrum at the iteration obtained from 

Eq. ([T|) using the recovered power spectrum P^^\k). Eq. ([3]), together with Eq. ([1]) for obtaining C;*-*'' from P^^\k) 
completely summarizes the standard RL method. The final recovered power spectrum is independent of the initially 
chosen P'^^^(fc) for reasonable forms of P^^)(fc) (see discussion in the Appendix of Ref. 

Due to noise and sample variance, the data Gf is measured within some non-zero error bars ct;. The standard RL 
method does not incorporate the error information at all and tends to iterate to fit features of the noise, as well. In 
our problem, the problem manifests itself as a non-smooth deconvolved spectrum P(fc) from the binned data that has 
poor likelihood with the full WMAP spectrum data. In Ref. [l| we devise a novel procedure to make the RL method 
sensitive to the errors (T; by modifying Eq. ([3]) to 



p^'+^\k) - p«(fc) = pW(fc) y G(/, k) 



G, 
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G, 



tanh 



(4) 



The idea is to employ a 'convergence' function to progressively weigh down the contribution to the correction p(*+i) — 
P''-* from a multipole bin where c|''* is close to GP within the error bar cr/. This innovation significantly improves 
the WMAP likelihood of the deconvolved spectrum. The final form of the recovered spectrum was obtained after 
smoothing the spectrum by a simple "bowler- hat" smoothing kernel, and varying width of smoothing of the kernel to 
get the best likelihood for the recovered angular power spectrum [33|. 

In Ref. we have carefully verified our improved Richardson-Lucy method on some toy models before applying 
it to the real data. The method has been shown to work well after a careful analysis of artifacts and convergence 
issues. Subsequently, a very similar primordial power spectrum has been obtained in Ref. [l9[ by deconvolution of 
the angular power spectrum from WMAP data using the method of regularized least squares. It is indeed of interest 
to apply other deconvolution methods and to compare the results. 
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FIG. 1: The recovered spectrum for a flat cosmological model with t—0.0, h—0.71, S7f,/i^ =0.0224 and f2A=0.73 (red-sohd 
Une). The green-dashed hne shows the reconstructed spectrum, using DWT, incorporating the coarse behavior and the low 
frequency features captured by the highest level wavelet coefficients. Even at this stage, the likelihood for the recovered angular 
power spectrum ( XeB = 974.1) is better than the likelihood for the best fit power-law spectrum ( XeS = 978.6). The blue- 
dotted line is the recovered spectrum using WMAP-3 data for the same background cosmology. This recovered spectrum gives 
AXe// = —15.93 with respect to the best power law primordial spectrum in the whole cosmological parameter space. 

III. DISCRETE WAVELET TRANSFORM 

Wavelet transforms provide a powerful tool for the analysis of transient and non-stationary data and is particularly 
useful in picking out characteristic local variations at different resolutions. This linear transform separates a data set 
in the form of low-pass or average coefficients, which reflect the average behavior of the data, and wavelet or high-pass 
coefficients at different levels, which capture the variations at corresponding resolutions. As compared to Fourier or 
window Fourier transform, wavelets allow optimal "time-frequency" localization simultaneously in the real, as well 
as, Fourier domain. The vocabulary of DWT stems from applications in one dimensional time-stream signal trains, 
but has found wide application in signal in other domains and dimensions. Specifically in our case, the 'signal' being 
transformed is the power spectrum, P{k), a one dimensional function of wavenumber, k. 

Wavelets are an orthonormal basis of small waves, with their variations primarily concentrated in a finite region, 
which makes them ideal for analyzing localized 'transient' signals. Wavelets can be continuous or discrete. In the 
latter case, the basis elements are strictly finite in size, enabling them to achieve localization, while disentangling 
characteristic variations at different frequencies [33 | . This is the primary reason for us to employ discrete wavelets for 
our analysis. 

In the construction of the basis set for discrete wavelet transform, one starts with the scaling function ip{x) (father 
wavelet) and the mother wavelet whose height and width are arbitrary: J ipdx = A, J ip dx = 0, J (ptpdx — 

0, / dx^l, J \tp\'^ dx = 1, where A is an arbitrary constant. In addition to the scaling and wavelet functions, their 
translates V'i.m = 2^/^^p{2^x — m), fj^m = 2^/'^ip{2^x — m), are also square integrable at different resolutions. Here, 
TO and j respectively are the translation and scaling parameters, both taking integral values with — oo < to < -l-oo. 
We start with the resolution value j — Q and increase it by integral units. The original mother wavelet corresponds to 
i/'o.Oi a-iid the father wavelet is given by (po,o- Higher values of j lead to the so called daughter wavelets which are of a 
similar form as the mother wavelet, except that their width and height differ by a factor of 2^ and 2^/^ respectively, 
at successive levels. The translation unit to/ 2^ is also commensurate with the thinner size of the daughter wavelet at 
resolution j . In the limit j oo, these basis functions form a complete orthonormal set. It needs to be pointed out 
that for a discrete data set with a finite number of points TV the maximum value of j is the largest integer less than 
or equal to logj N. 
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A signal f{x) can then be expanded as 

+00 +00 

f{^)= X! '^O.mVQ,m{x) + ^ ^dj^rn-lpj^rn{x) . (5) 
m— — 00 rn— — ooj>0 

Here, c^^m's are the low-pass coefBcients and dj.m's are the high-pass or wavelet coefficients, respectively capturing 
the average part and variations of the signal at resolution j and location m. In practice, for a finite data set one starts 
with the highest level of resolution and progressively moves to resolution on grosser scales keeping in mind the physics 
of the problem and the size of the data set. The lowest level of resolution one chooses then corresponds to j = and 
the higher levels correspond to larger values of j. For the discrete wavelets, the property of multi-resolution analysis 
(MRA) leads to Cj^m = X^n '^('^ ~ 2m)cj+i^„ and dj ^ = J2n ^("' ~ 2TO)cj_|_i_„, where h{n) and h{n) are respectively 
the low-pass (scaling function) and high-pass (wavelet) filter coefficients, which differ for different wavelets. Thus, 
both low-pass and high-pass coefficients at a resolution j can be obtained from the low-pass coefficients at a higher 
resolution j + 1- The low-pass coefficients Cj+i^m are obtained though the convolution of the signal f{x) with the 
scaling function (fj+i^m = 2(-'+^^/^<y9(2-'+^a; — to). For a fixed m, in the limit j oo, the scaling function becomes a 
Dirac delta function and hence the corresponding low-pass coefficient is the signal itself at point to. This implies that, 
starting from the finest resolution of the signal, one can construct both scaling and wavelet functions, by convolution 
with the corresponding filter coefficients. Hence one can carry out the wavelet decomposition, as also the inverse 
transform, with the help of h{n) and h(n), without explicitly knowing the basis set. In this respect, wavelet transform 
is significantly different from the Fourier transform. 

For the Haar wavelet, h{Q) = h(l) = and h{Q) — ~h{l) = The Haar basis is unique, since it is the only 
wavelet which is symmetric and compactly supported. In a level one Haar wavelet decomposition, the level-I low-pass 
(average) and high-pass (wavelet or detail) coefficients are respectively given by the nearest neighbor averages and 
differences, with the normalization factor of In the subsequent step, the average coefficients are divided into 
two parts, containing level-II high-pass and level-II low-pass coefficients. The high-pass coefficients now represent 
differences of averaged data points corresponding to a window size of two. Thus higher level coefficients represent 
lower frequency features. Wavelets belonging to the Daubechies family are designed such that the wavelet coefficients 
are independent of polynomial trends in the data. We have carried out a 10-level decomposition using Daubechies-4 
wavelets for isolating fluctuations at different resolutions. Daubechies-4 wavelet satisfies / xip{x)dx — 0, in addition to 
all other conditions. Because of this the wavelet coefficients capture fluctuations over and above the linear variations. 
Simple models of inflation predict that log P{k) vs log(fc) will have a linear relation. The choice of Daubechies-4 
wavelets is the minimal (simplest) wavelet that allows us to study variations about linear behavior in a window whose 
size increases with the level of decomposition. 

Furthermore, to study the significance of features located at different wavenumbers and at different resolutions, 
fluctuations associated with wavelet coefficients of different levels are added to the average behavior captured by the 
low-pass coefficients in order to reconstruct a smoothened power spectrum. Then a likelihood analysis with respect 
to the WMAP data is performed. 

If the data set is of size 2^ a maximum of N level decompositions can be carried out. In the case of fewer data 
points, one needs to supplement the data with additional points to carry out an N level decomposition. Due to the 
finite size of the filter coefficients, one also encounters boundary artifacts due to circular or other forms of extensions. 
In our case, for minimizing these boundary artifacts, we have carefully padded the data with constants at both ends, 
which were then removed after reconstruction. We worked with a data set of 8192 points for log P(k) vs log(/c) 
recovered from WMAP-1 observations. The data set was arranged to be equally spaced, as required by wavelet 
programs. 

IV. FEATURES OF THE PRIMORDIAL POWER SPECTRUM 

One of the most challenging questions of the modern cosmology is to find an inflationary scenario satisfying all the 
cosmological observations. The shape of the primordial power spectrum has the key role in this investigation. In our 
previous paper 1] we have used Richardson-Lucy deconvolution algorithm to find the shape of the primordial power 
spectrum using the cosmic microwave background data. 

By using DWT, we decompose the recovered primordial power spectrum locally at different resolutions. We then 
calculate the angular power spectrum after including the variations in the primordial power spectrum at different 
resolutions. We subsequently compute the likelihood of the primordial power spectrum at each of these stages so as 
to quantify the effect of different features on the recovered angular power spectrum. We use WMAP-1 likelihood code 
available at the LAMBDA site [1^ and quote our likelihood C in terms of xle = — 21n£. 
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FIG. 2: The left panels show the DWT decomposition of the features of the primordial power spectrum. Step A shows the 
reconstructed primordial power spectrum using the 10th level low-pass coefficients. Steps B, C, D, E and F show the localized 
variations due to the wavelet coefficients at the 10th, 9th, 8th, 7th and the 6th levels respectively. These variations in the 
primordial power spectrum are most prominent close to the horizon scale and are significant only in the first few panels, 
corresponding to 'low frequency' variations. The panels on the right compare the resultant angular power spectrum, CI with 
the binned WMAP-1 angular power spectrum data. Going down from the top we progressively add features from different 
levels of wavelet coefficients (left) to the primordial power spectrum and show the difference, Ci — Cf and error bars for Cf 
(right). 



The original primordial power spectrum which we use in this paper (see figure [T|) is the final recovered primordial 
power spectrum obtained in Ref. [l| for a flat ACDM cosmological model with r=0.0, /i=0.71, rif,/i^=0.0224 and 
nA=0.73. The resultant Ci spectrum using this recovered spectrum has a likelihood far better than a scale invariant 
or a best fit power law spectrum. In this paper, we investigate how the features of this spectrum contribute to 
improving the likelihood. 

The blue dotted line in figure [1] represents the recovered primordial spectrum by using WMAP-3 data which gives 
Axgj^ = —15.93 (by using WMAP-3 likelihood code) with respect to the best power law primordial spectrum in 
the whole cosmological parameter space. In this paper we are not using this recovered spectrum for our wavelet 
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FIG. 3: The primordial power spectrum reconstructed with the low-pass coefficients is shown in the top panel. Addition of 
features of steps B, C and D is shown in the middle panel. The combination of both these panels is shown in the bottom 
panel. Note the significant effect of the features on the hkelihood. 



analysis and it has been presented here just to show the close similarity of the features of the recovered spectrum 
from WMAP-1 and WMAP-3 data. 

First, we smooth the spectrum using DWT, and then we systematically include features at different resolutions and 
calculate the likelihood with respect to the WMAP data. 

In figure [21 left panel (step A), we show the coarsest behavior of the data reconstructed using only the 10th level 
low-pass coefficients. The right panel shows the resultant compared with the observed binned data Cf, and its 
error bars. The likelihood of the at this stage corresponds to Xcs ~ 994.4 which is better than the Xcfr for the 
H-Z spectrum. In the right panel (step B) we see that if we include 'low frequency' features as captured by the 10th 
level high-pass coefficients shown in the left panel (step B) to the previous spectrum at step A, the resultant Cf will 
be closer to the observed data. Hence we expect the likelihood to be improved. In fact this improves the likelihood 
significantly and the xls = 974.1 at this stage (the reconstructed spectrum is the green-dashed fine in figure [T]). This 
xlff is better than the best fit power law spectrum with the xte = 978.6. 

As we progressively add back more features at 'higher frequencies' to the spectrum the likelihood improves signifi- 
cantly. The 9th, 8th, 7th and 6th level wavelet coefficients are shown in the left column of figure^] (steps C, D, E and 
F) . The plots in the right panel of figure [2] show how we get closer to the observed data within the error bars as we 
include more wavelet coefficients. The likelihoods at these steps correspond to are x^^^ = 964.8 at step C, = 961.0 
at step D, Xcff = 957.4 at step E, and Xcfr = 956.5 at step F. 

In figure [3] we compare the coarsest spectrum of step A and the spectrum after adding the features of step B, C 
and D. The local features which are responsible for the significant improvement of the likelihood are clearly seen. 
We note that most of these features are localized around the horizon scale. In figured! we see that after few stages 
of adding the more detailed features to the spectrum, the likelihood does not improve anymore and agreement with 
the observed data is not sensitive to these features. 



V. CONCLUSION 



This paper presents a detailed analysis of the recovered primordial power spectrum for a fiat ACDM cosmological 
model, given in Ref. UJ. The recovered spectrum has a likelihood far better than a simple scale invariant or a scale 
free spectra. In Ref. [I^l, similar features for the primordial spectrum have been detected by using a completely 
different method and the significance of the features have also been evaluated. In this paper we use Discrete Wavelet 
Transform to decompose the features of the spectrum to quantify and understand their role in improving the likelihood. 
In addition to the infra-red cut-off around the horizon, which was proposed by many groups to explain the lack of 
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FIG. 4: The likelihood values improve as detailed features are added to the primordial power spectrum reconstructed with 
low-pass coefficients. A corresponds to the primordial power spectrum reconstructed with low-pass coefficients only. B — K 
correspond to the primordial power spectrum with the contribution of wavelet coefficients of levels 10-1 progressively added to 
the coarse spectrum A. There is no significant improvement in the likelihood beyond level 6 (step F). Green-dashed lines mark 
the Gaussian equivalent a-levels of the likelihood relative to the best recovered spectrum. Results are compared with the best 
likelihood given by Power Law and Harrison-Zeldovich spectrum. Nc is the number of coefficients (low-pass and high-pass) 
used to define the primordial spectrum at each step. 



power in very low multipoles of the observed angular power spectrum, we show that the features around the horizon 
are playing a crucial role in improving the likelihood. In fact, the effect of these features on improving the likelihood, 
are very significant (figure [3]). We find that these strong features are localized around the horizon. Work on WMAP 3 
years data is in progress and will be reported in a forthcoming publication [181] . The work in progress implies refinement 
of methodology includes directly employing DWT right at the stage of smoothing the raw recovered spectrum and 
dispensing with the need to remove a 'template' spectra as in Ref. [ij. 
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